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Abstract 

In this note we provide explicit expressions and expansions for a special function J which 
appears in nonparametric estimation of log-densities. This function returns the integral of 
a log-linear function on a simplex of arbitrary dimension. In particular it is used in the R- 
package LogCondDEAD by Cule et al. (2007). 
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1 Introduction 

Suppose one wants to estimate a probability density / on a certain compact region C C M"', based 
on an empirical distribution P of a sample from /. One possibility is to embed C into a union 

m 

s =[js, 

of simplices Sj C M.'^ with pairwise disjoint interior. By a simplex in we mean the convex 
hull of d + 1 points. Then we consider the family Q = Q{Si , ■ ■ ■ , Sm) of all continuous functions 
: 5 — > R which are linear on each simplex Sj. Now 

ip := argmaxl / ipdP— / exp(^/'(x)) I (1) 
i'eg \Js Js J 

defines a maximum likelihood estimator / := exp('(/;) of a probability density on S, based on P. 
For existence and uniqueness of this estimator see, for instance, Cule et al. (2008). 

To compute t/j explicitly, note that ip ^ G is uniquely determined by its values at the corners 
(extremal points) of all simplices Sj, and f ip dP is a linear function of these values. The second 
integral in ^ may be represented as follows: Let Sj be the convex hull of XQj ,xij, . . . , Xdj S W^, 
and set yij := 'ilj{xij). Then 



„ m „ m 

j exp(^/;(x))dx = exp(V'(x)) = ■ J{yoj,yij, 
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where 

Dj := det[xij - xoj,X2j -xoj,... , xaj - xoj] , 
while J(-) is an auxiliary function defined and analyzed subsequently. 



2 The special function J{-) 
2.1 Definition of J( ) 

For d G N and r > let 

Td{r) := {ue[0,^y ■.Y,u^<ry 
Then for yQ,yi, ... ,yd E M we define 
Jiyo,yi, ■ ■ ■ iVd) ■= 



1=1 



/ exp({l -u+)yo + y^Uiyi]du 
Mil) V ^ y 



with u+ := Y^i^i Ui. 

Standard considerations in connection with beta- and gamma-distributions as described in Sec- 
tion |6]reveal the following alternative representation: 

1 

Jiyo,yi,---,yd) ■= ^Eexp(^^5iyi^ 

i=0 

with Bi = Bfi.i := Eij Ylt=o-Es and stochastically independent, standard exponential random 
variables Eq, Ei, . . . , E^. This representation shows clearly that J(-) is symmetric in its argu- 
ments. 

An often useful identity is 

J{yo,yi,---,yd) = exp{y^)J{yo-y^,yi-y^,...,yd-y*) for any G M. (2) 

2.2 Some useful recursions 

It is well-known that for any integer < j < d, 



Ei \ ^ _ YJi=QEi ( E, 



V Ei=o Es J Es=o Es V Tls=j+i Es J .^.^^ 

are stochastically independent with B Beta(j + 1, d — j); see also Section[6] Hence we end up 
with the following recursive identity: 



Jiyo,yi, ■■■,yd) 
jKd-j-iy. 



E( J(i?yo, • • • , Byj)J{{l - B)yj+i, . . . , (1 - B)yd)) 







1 

u^{l - u)'^'^'^J{uyo, . . .,uyj)J{{l - u)yj+i, . . . , (1 - u)yd)du 
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with 



J(r) := exp(r). 



Here we used the well-known identity 

/(I - ufu"^ du = 



£lm\ 



for integers i,m > 0. 



{e + m + iy. 

Plugging mj = d—l into the previous recursive equation leads to 

Jiyo,yi,---,yd) = / u'^~^J{uyo,...,uyd-i)exp{{l-u)yd)du. 
Jo 

For d = I these recursive formulae are useless, but one compute J(yo, yi) explicitly: 



(3) 



(4) 



J{yo,yi) 



exp 



((1 - u)yo + uyi) du 



( exp(yi) - exp(yo) , 

- II yo 7^ y\ 



y\ - yo 

exp(yo) 



if yo = yi- 



But now we can prove an alternative recursive formula: For arbitrary integers d > I and real 
numbers yo, yi, • • • , yd, 



J{yo,yi, ■■■,yd) 



J{yi,y2, ■■■,yd)- J{yo,y2, ■■■,yd) 



d 



if yo / yi, 
yi - yo 



■^J{yi,y2, ■■■,yd) 
dyi 



if yo = yi- 



(For d = I, just erase "2/2, • • •" on the right hand side.) We prove this formula by induction on 
d and restrict our attention to the case yo 7^ yi- The case yo = yi follows by a suitable limit 
argument. Suppose that ^ is true for some d > 1. Then it follows from the induction 
hypothesis and a second application of & that 

<^(yo,yi, • • ■,yd,yd+i) 
1 

u'^Jiuyo,uyi,-- ■ ,nyrf)exp((l - u)yd+i) du 

^ dJ(.uyi,uy2,...,uyd)-J{uyo,uy2,...,uyd) . 

u" exp((l - u)yd+i) du 

uyi - uyo 

^ d~i J{uyi,uy2,.. ■,uyd) - J{uyo,uy2,.. .,uyd) 



yi - yo 



exp((l - u)yd+i) du 



yi 



(/ ^J{uyi,uy2,...,uyd)exp{{l-u)yd+i)du 

- yo ^Jo 



u"^ V(Myo,My2, . . . ,nyd)exp((l -u)yd+i)du 



Jiyi,y2, ■ ■ ■,yd,yd+i) - J{yo,y2, ■ ■ ■,yd,yd 



+1; 



yi - yo 
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3 An expansion for J(-) 

With y := {d+ Yli=o Vi ^'^'^ '■= Vi - V one write 

J{yo,yi,---,yd) = ex.p{y)J{zo,zi,...,Zd) 
by virtue of Note that z+ := Ya=o = 0- Thus, as z := {zi)f^Q 0, 

d 1 

d\ J{zo, zi,...,Zd) = l + Yl + 2 5Z HBiBj)ziZj 

i=0 i,j=0 

1 

+ -^Y1 ^{B,B^Bk)ziZjZk + 0{\\z\ 

It follows from Lemma [6?T] that 
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i,j,k=0 



^i^) = n ^^'/t*^ + ^niegers ko,ki,...,kd>0. 

i=0 1=0 



In particular, 



d+l 



IE(S3) = T— TEiB^Bi) = jj^, JE{BoBiB2) 



Consequently, Eto^C^O^i = ^{Bq)z+ = 0, 

[d + 2]2 ^ W.{B,Bj)ziZj = (l{i = j} • 2 + l{i / 

= 5Z(l{i=j} + l)2i^, 



d 



i=0 
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and 



[d + 3]3 ^{BiBjBk)ziZjZk 

i,j,k=0 
d 

= {ni=j = k]-Q + l{Mh3,k] = 2]-2 + l{#{i,j,k]=i})ziZjZk 

i,j,k=0 
d 

= Y =j = k}-5 + = 2} + l)z,zjzk 

i,j,k=0 

d d 

= 5^zf + 3 ^ l{s^t}z^^zt + zl 

i=0 s,t=0 
d d d 



i=0 
d 



i=0 s=0 s=0 

i=0 



Consequently, 

i=0 i=0 

4 A recursive implementation of J(') and its partial derivatives 

By means of ^ and the Taylor expansion Q one can implement the function J(-) in a recursive 
fashion. In what follows we use the abbreviation 



Va-A 



iya,---,yb) ifa<6 
if a > 6 



To compute J{yo:d) we assume without loss of generality that yo < yi < ■ ■ ■ < Ud- ^t. follows 
from (|5]l and symmetry of J(-) that 

J. . J (yi-.d) - J (yo-.d-i) 

J[yo:d) = 

yd - yo 

if yo 7^ yd.. This formula is okay numerically if y^ — yo is not too small. Otherwise one should use 
This leads to the the pseudo code in Table [T] 

To avoid messy formulae, one can express partial derivatives of J(-) in terms of higher order 
versions of J(-) by means of the recursion For instance, 

dJ{yo:d) _ J{yo + e, yy.d) - Jjyo, yi-.d) 



dyo e^o 
= lim 

= Jiyo,yo,yi: 



lim J{yo,yo + e, yi-.d) 
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Algorithm J ^ J{y,d,e) 
if yd - yo < e then 

^2-Eto(2/^-y)V2 

^3-Eto(y-y)V3 

J ^ exp(y)(l/d! + Z2/{d + 2)! + Z3/{d + 3)!) 
^ {3{yi:d,d- l,e) - J(yo:d-i,d - l,e))/(yrf - Vo) 



else 



end if. 



Table 1: Pseudo-code for J{y) with ordered input vector y. 



Similarly, 



lim( 



J {yo + e, yi:d) - J {yo, yi-.d) J {yo, yi-.d) - J{yo- e, yi-.d) 



2 lim 



J {yo, yo + e, yi:d) - Jjyo^yo- e, yi:d) 
2e 



2 lim J{yo,yo -e,yo + e, yi:^) 

e— >0 

2 J(?/o,yo,yo, 



while 



dyodyi 



lim f 



(j/o + e, yi + y2:d) - (j/o, yi + e, y2:d) 
e 



^.^ Jjyo^yo + g,yi + e, y2:d) - J{yo,yo + ^2:0 
lim J{yo,yo + e,yi,yi + e, y2:d) 

e— »0 

J{yo,yo,yi,yi, y2:d)- 



5 The special cases d = 1 and d = 2 

For small dimension d it may be worthwhile to work with non-recursive implementations of the 
function J(-). Here we collect and extend some results of Diimbgen et al. (2007). 



5.1 General considerations about a bivariate function 

^ M which is infinitley often differentiable. 

f{s) - fir) 



In view of ^ we consider an arbitrary function / : 
Then 



h{r,s) := 



s — r 



fir) 



fir) + ^is-r)+0{{s-r)^) 



if s 7^ r, 



as s — > r, 
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defines a smooth and symmetric function h : M? 
two are given by 



dh{r, s) 
dr 



d'^h{r, s) 

Qj.2 



drds 



R. Its first partial derivatives of order one and 

f{s)-f{r)-f{r){s-r) 



{s - ry 



if s 7^ r, 



Z!|:)+Z!^(,_,) + 0((,_,)2) as.^r, 
f 2{fis)-f{r)-f'{r)is-r))-is-rff"ir) . 



= < 



(s — r)^ 



if s 7^ r, 



f"'ir) , f""ir) 



H T7^(* - r) + 0((s - r) ) as s ^ r, 



3 12 

r (5-r)(/'(r) + /'(^))-2(/(^)-/(r)) 



(s — r)^ 



if s ^ r, 



f"{r) ^ _ ^) + o ((^ _ ^)2^ ^ 



r. 
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The other partial derivatives of order one and two follow via symmetry considerations. 

5.2 More details for the case d = 1 

Here the general formula boils down to 



J{r,s) = / exp((l — u)r + us) du, 
Jo 



and elementary calculations show that 



J(r,s) = 



exp(g) - exp(r) .^^^^^ 
exp(r) if r = s. 



This is just the function introduced by Diimbgen, Hiisler and Rufibach (2007). Let us recall some 
properties and formulae for the corresponding partial derivatives 

Qa+b rl 

Ja,b{r,s) := Q^aQgb ^) = {l-uyu^exp({l-u)r + us)du. 
Note first that 

Ja,b{r,s) = Jb,a{s,r) = exp{r)Ja,b{0,s -r). 
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Thus it suffices to derive formulae for (r, s) = (0, y) and b < a. It follows from ([3]l tliat 



In particular, 



-^,0(0,^) 



•^'J fc=0 

(l-n)Vdn./ 

fc=o ■ 

00 1 
\ - a! 

-^(exp(y)-^l^ 



-/i,o(0,2/) 



^2,0(0, y) 



^3,0(0,2/) 



J4,o(0,2/) 



exp(y) - 1 - y 

2(exp(7/) - 1 - - yV2) 

1 + A + y! + ^ + 0(/) (y^O), 
3 12 60 360 ^ ^ ^ ^ 

6(exp(y) -l-y- yV2 - j/V6) 
24(exp(j/) - 1 - 7/ - 2/V2 - yV6 - yV24) 



5 30 210 1680 
Another general observation is that 

Ja,b{r,s) = [ {l-uy{l-{l-u))''exp{{l-u)r + us)du 
Jo 



i=0 

In particular, 



5^(^)(-l)V.+,,o(r,.). 

i—n V / 



Ja,iir,s) = Ja,oir,s) - Ja+iflir,s), 

Ja,2{r,s) = Jafi{r,s) -2Ja+i,Q{r,s) + Ja+2fi{r,s). 
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On the other hand, 



k=o ■ 

oo 

= 



k=0 



a\[k + b]b ^ 
{k + a + h+l)\^ 



with [r]o := 1 and [r\m '■= Y\a!=o^{''^ ~ for integers m > 0. In particular, 



Ji,i(0,y) 

5.3 The case d = 2 

Our recursion formula ([5]) yields 

J(r,s,t) 



exp(y)(y-2) + 2 + y 



J(.,t)- J(r,t) 
II r 7^ s, 

s — r 

Jio(r,t) ifr = s. 



Because of J's symmetry we may rewrite this in terms of the order statistics y(o) < < 2/(2) of 
(yi)Lo as 

' J{y(i),y(2)) - ^(y(o),2/(i)) 



Jir,s,t) 



y{2) - y(o) 

exp(y(o)) 



if 2/(0) < y(2), 



if 



y(o) = 2/(2)- 



For fixed third argument t, this function J(r, s, t) corresponds to h{r, s) in Section [5?T] with 
/(x) := J(x, t). Thus 

J(s,t) - J(r,t) - Ji,o(r,t)(s-r) 



dJ{r, s, t) 
dr 



(s — r)2 
J2fl{r,t) J3fi{r,t){s -r) 



+ 
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if r / s, 



+ 0((s-r)^) as s ^ r. 



Moreover, 

d^J{r,s,t) 

d^J{r,s,t) 
drds 



2{j{s, t) - J{r, t) - Jifl{r, t){s - r)) - {s - rfJ2,o 



(s — r)3 
'/3,o(^0 J4,o{r, t){s - r) 
3 12 

( Ji,o(r, t) + Ji,o(s, t)) (s - r) - 2( J(s, t) - J(r, t)) 



if r 7^ s, 



(s — r)^ 



+ 0((s — r)^) as s ^ r, 

if r 7^ s, 



6 



+ 



12 



+ 0((s-r)2) 



as s — > r. 
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6 Gamma and multivariate beta (Dirichlet) distributions 

Let Go, Gi, . . . , Gm be stochastically independent random variables with Gi ~ Gamma(ai) for 
certain parameters Oj > 0. That means, for any Borel set A C (0, oo), 

P(G, e^) = / riai)-'y''^-^expi-y)dy. 
J A 

Now we define a^. := X^I^o ^+ X^ilLo ^'^'^ 

-B := (Gi/G+)™0' ^ '■= iGi/G+)'^i. 
Note that B is contained in the unit simplex in M"*+^, while B is contained in the set 

%n := {ue{0,ir -.u+Kl} 
with u+ := YliLi ^i- We also define uq := 1 — u+ for any u e 7^. 

Lemma 6.1 The random vector B and the random variable G+ are stochastically independent. 
Moreover, 

G+ ^ Gamma(a+) 
while B is distributed according to the Lebesgue density 

on 7^. For arbitrary numbers kQ, ki, . . . , km > and k^ := X^i^o 



IE 



rn 

(n 

i=0 



r(a+) 



r(a+ + k+) jLl r(ai) 



n 



As a by-product of this lemma we obtain the following formula: 
Corollary 6.2 For arbitrary numbers oq, ai, . . . , > 0, 

/ n<^^'^" = r(a+)-inrK) 



1=0 



Proof of Lemma 16^1 Note that G = {Gi)'^Q my be written as H(G+,-B) with the bijective 
mapping H : (0, cxo) x 7^ — > (0, oo)"'^^, 



Note also that 



det DE{s,u) = det 



/ Uq —s —s 

ui s 

U2 s 

\um • • • 



: 
s J 



det 



/ 1 

ui s 

U2 S 



0\ 





\um ••• sy 



11 



Thus the distribution of B) has a Lebesgue density h on (0, oo) x %n which is given by 

m 

h{s,u) = JJ(r(ai)"-^H(s,'u)°*""^exp(-H(s,it)i)) • |det 1)5(5, 

i=0 
m 

= l[{T{a,)-\suir-' expi-sui)) ■ s"^ 

1=0 

m 

= T{a+)-h-+-'exp{-s)-f{u). 

Since this is the density of Gamma(a+) at s times f{u), we see that G+ and B are stochastically 
independent, where G+ has distribution Gamma(a+), and that / is indeed a probability density 
on Tm describing the distribution of B. 

The fact that / integrates to one over 7^ entails Corollary 16.21 But then we can conclude that 

TTi Tfi TTl 

i=0 "^^"i i=0 ^™ i=0 

^ r(a+) -pr r(aj + kj) ^ 
r(a+ + A;+)ii r(a,) ' 
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